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ABSTRACT 


Martin et al. 


(2014b) showed that a substantially misaligned accretion disk 


around one component of a binary system can undergo global damped Kozai- 
Lidov (KL) oscillations. During these oscillations, the inclination and eccentricity 
of the disk are periodically exchanged. However, the robustness of this mechanism 
and its dependence on the system parameters were unexplored. In this paper, 
we use three-dimensional hydrodynamical simulations to analyze how various 
binary and disk parameters affect the KL mechanism in hydrodynamical disks. 
The simulations include the effect of gas pressure and viscosity, but ignore the 
effects of disk self-gravity. We describe results for different numerical resolutions, 
binary mass ratios and orbital eccentricities, initial disk sizes, initial disk surface 
density profiles, disk sound speeds, and disk viscosities. We show that the KL 
mechanism can operate for a wide range of binary-disk parameters. We discuss 
the applications of our results to astrophysical disks in various accreting systems. 


Subject headings: accretion, accretion disks - binaries: general - hydrodynamics - 
planets and satellites: formation 
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Introduction 


Disks in binary systems may be misaligned with respect to their binary orbit planes. 
Alignment is expected in cases where the disk formation process itself is coplanar, such as 
the case of mass-exchange binaries. In that case, a coplanar gas stream originates from 
the inner Lagrange point of the mass-losing star and results in the formati on of a cop la nar 


disk. Coplanarity is expe cted if the binary forms from disk fragmentation (IBonnell fe Bate 


misa 

igned disk, c 

ue to 

the effects of tidal dissir 

)ation in the disk 

(Papaloizou & Terauem 

1995 

Bate et al. 

2000; 

Lubow & Ogilvie 

2000; 

King et al. 

2013 

■ 


On the other hand, disk misalignment is expected to persist for some time period in 
cases where the binary forms from a merger pro cess , such as occurs with supermassive black 


hole (SMBH) binary systems (e.g., 


King & Pringle 


2006). If a young binary star system 


accretes material after its formation process, the material is likely to be misaligned wi 


binary orbit and so misaligned disks may form around young stars (Bate et ah 


h the 


2010 ). In 


addition, a co planar dis 


warping (Pringle 


1996 


c may become noncoplanar due to an instability , such as radiation 


Wijers & Pringle 


1999 


Ogilvie & Dubus 


2001). Be/X-ray binaries 


may_ ha ve t ransient disks that are misaligned with respect to the binary orbit flMartin et ah 


2009 


2011). These disks may be produced by the expulsion of equatorial material from a 


rapidly r otati ng Be star whose spin axis is misaligned with the binary rotation axis (e.g., 


Porter fe Rivinus 


200 3). 


There is both direct and indirect observational evidence for disk noncoplanarity. For 
widely separated stars in a binary, greater than 40 AU, a misaligned disk may occur 
because the stellar equatorial inclinations, based on spin s, are ob servationally inferred to be 


misaligned with respect to the binary orbital planes (Hale 


1994 ). The young binary system 


HK Tau provides direct evidence for noncoplanarity, since disks are observed around both 
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components with one disk edg e-on and the other more face-on (jStapelfeldt et al. Il998h . 


Recent ALMA observations by 


Jensen fe Akeson 


(120141) suggest that its disks are mutually 


misaligned by more than 60°, although the plane of the binary orbit is not known. Strong 


mutual disk misalignment (~ 72°) has lately a 


so been detected for the two circumstellar 


disks in V2434 Ori. a binary svstem in Orion (Williams et ah 

2014). 

Less direct evidence 

comes from the existence of extra-solar planets w 

lose orbits are tiltec 

with respect to the 

soin axis of the central star (Albrecht et al. 2012; 

Winn & Fabrvckv 

2014 

). If the planets 

reside in binary star systems, this evidence suggests that these planets may have formed in 


disks that are misaligned with the binary orbital plane (e.g., 


Bate et al. 


2010 


2012 ). 


Batygin 


Test particles whose orbits are sufficiently inclined with re spect to the_ 


circular orbit binary can undergo Kozai-Lidov (KL) oscillations (jKoza i 


1962 


plane o 


Lidov 


1962) 


In these oscillations, the particle’s orbital inclination and eccentricity evolve in such a 
manner that inclination is exchanged for eccentricity. Due to this process, a test particle 
orbit that is initially circular can achieve a high eccentricity when it evolves to smaller 
inclination. For example, an initially circular orbit at an inclination of 60° achieves an 


effects can occur for eccentric orbit binaries ( 

7 ord e 

al. 

2000 

; Lit 

hwick & Naoz 

2011 

Naoz et al. 

2011, 

2013a 

b; 

Tevssandier et al. 

2013; 

Li et ah 

2014 

Liu et al. 

201 

5). 


The KL effect for ballistic objects has been applied to a wide range of astr onomical 


pr ocesses. These inclu de inclinations of asteroids and irregular satellites (IKozai 


1962 


Nesvornv et al. 2003) 

big 

ers ( 

i orbital eccentricities of some extra-soh 

rr planets and 

formation of hot Jupil 

dolman et al. 

1997 

Mazeh e 

al. 

1997; 

Wu & Murray 


2003; 

Takeda & Rasio 

2005 


Xiang-Gruess & Papaloizou 

2013; 

Dawson & Chians: 

2014; 

Done et ah 20lJ: Petrovich 

2015; Rice 

2015) 

misalignment between an exoplanet’s 
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orbital axis and stellar spin axis 

(Storch et ah 

2014, 2015), black hole mergers (Blaes et al. 

2002: 

Miller & Hamilton 

2002) 

, the formation of close binary stars 

(Harrington 1968; 

Mazeh & Shaham 

1979; 

viscleva et al. 1998: 

Fabrvckv & Tremaine 

2007), tidal disruption 

events (Chen et al. 

2011 

), and the formation of Type la supernovae 

Kushnir et al. 2013). 


In a recent Letter ([Mart in et al. 


2014b . hereafter Paper I), we showed in our Smoothed 


Particle Hydrodynamics (SPH) simulations that a fluid disk with pressure and viscosity 
that orbits a member of a binary can also undergo global KL oscillations. In this paper, 
we explore this process in more detail by considering how it is affected by various disk and 
binary parameters. Following the approach of Paper I, we attempt to relate the properties 
of disk KL oscillations to the properties of KL particle oscillations. 


The outline of the paper is as follows. In Section [2] we describe the properties of 
particle orbits in binaries. The parameters of the particle orbits and the binaries are chosen 
to be similar to those of the disk simulations. In Section [3] we describe the results of 
SPH simulations for a range of initial particle numbers, disk viscosities, disk aspect ratios, 
initial density profiles (including changes to initial inner and outer disk radii), initial disk 
inclinations, binary mass ratios, and binary eccentricities. Section [4] contains the discussion 
and summary. 


2. KL cycle of a test particle 

As discussed in Paper I, the properties of a fluid disk undergoing KL oscillations are 
related to the properties of test (massless) particles undergoing KL oscillations that orbit 
at similar radii. In both the disk and test particle cases, the inclination and eccentricity are 
periodically exchanged. However, there are some fundamental differences. The period of a 
test particle oscillation varies with distance from the central object. In the case of a disk 
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that is subject to pressure and viscous forces, the oscillations are global (i.e., the period 
is the same over all radii) and subject to damping by dissipation, unlike the test particle 
case. The global disk oscillation frequency is expected to involve an angular momentum 
weighted average of the test particle oscillation frequencies (Equation (4) of Paper I). 
Before considering the disk oscillations, we describe in this section the behavior of test 
particle orbits. We apply binary parameters that are similar to those we adopt for the disk 
simulations described in Section [3j 

We consider ballistic particle orbits around one component of a binary system. Each 
particle orbit is initially circular, but substantially inclined with respect to the orbital 
plane of the binary. If the binary orbit is circular, then over long timescales, the binary 
companion effectively acts as a uniform circular ring. Since this effective potential is time 
independent, the energy of the particle about the central mass is conserved. Thus, the 
semi-major axis a of the particle in a binary is nearly independent of time t, that is 

a(t) ~ constant, (1) 

where the constant is determined by initial conditions. Because the time-averaged 
perturbing potential is axisymmetric, the component of the angular momentum of the 
particle that is perpendicular to the binary orbital plane is approximately conserved over 
long timescales. From this conservation principle and Equation (JTj) we have that 

\/l — e 2 (t) cos (i(t)) ~ constant, (2) 

where the constant is determined by initial conditions, e is the orbital eccentricity of the 
particle, and i is the inclination of the particle orbit with respect to the binary orbital plane. 

We denote the masses of the binary component objects as M c and M p , where M c (M p ) 
is the mass of the central (perturbing) object to the particle orbit. Analytic calculations 
have determined that the KL oscillation period for a particle subject to the effect of a 



binary in an eccentric orbit is approximately given by 


t K l 

P b 


M r 


M p P h 


Mr, 


P 


(1 - 


( 3 ) 


(Holman et al. 


1997 


Trimmer] et ah 


1997|; 


Kiseleva et al. 


1998 b where P = y/GM c /a 3 is 


the orbital period of the particle with orbital semimajor axis a, Pb = 2tt /fib is the orbital 
period of the binary, eb is the eccentricity of the binary orbit, and fib = \JG(M C + M p )/a\ 
is the binary orbital frequency for binary semi-major axis Ob- 


For the KL oscillations to occur, the initial inclination of the test particle orbit Zo 
must satisfy the condition that cos 2 (z 0 ) < cos 2 (z cr ) = 3/5, where z cr is the critical angle 
for KL oscillations. This condition requires that 39° < ?o < 141°. During the KL cycle, 
the inclination i of an initially circular prograde particle orbit (zo < 90°) oscillates between 
the initial value z 0 and the critical angle z cr ~ 39°, while the eccentricity oscillates between 
e m in = 0 and the maximum value of 


6 max 



5 

3 


cos 2 (z 0 ), 


( 4 ) 


that is achieved for i = z cr , w' 
constant value of cos (z 0 ) (e.g., 


ri ch follows from Equation ([2]) with the right-hand side 


Trimmer] et al. 


1997) 


For an eccentric binary orbit, the vertical component of a particle’s angular momentum 
is not conserved. Equation (J2]) breaks down o ver long timescales and the maximum 


eccentricity of the particl e can approac 

”3 


Naoz et ah 


2011 


2013al bl: 


Liu et. al. 


l unity (Ford et al. 


2000 


Lithwick fe Naoz 


2011 


2015). The orbit can flip its orientation from 


prograde to retrograde with respect to the binary. This flip occurs as the orbit passes 
through a radial state (e = 1). In passing through a nearly radial orbit, the particle may 
collide with the central object. 


In our numerical calculations, we compute the eccentricity of a test particle with 
coordinates r(t) by considering its specific angular momentum j(t) relative to the disk 















































central object at position r c (t), 


j(t) = (r — r c ) x (r — r c ), 


(5) 


where the dot denotes differentiation in time. The relative specihc energy of the particle is 


E(t) = -\r - r c 


GM C 


\r — r r 


and thus its eccentricity is 


e (t) ~ 1/1 + 


mj? 

0 GM c y 


The evolving inclination of the particle is 


i(t) = arccos [ ) , 

\3\ 


( 6 ) 


(7) 


( 8 ) 


where j z is the component of the angular momentum in Equation ([5]) that is perpendicular 
to the orbital plane of the binary. 


Figure |Tj plots the numerical results for a test particle undergoing KL oscillations 
from simple three-body simulations. The upper and lower rows show the evolution of the 
particle’s orbital eccentricity and inclination, respectively. We investigate the effects of 
changes to four system parameters: the initial orbit inclination i 0 , the initial orbital radius 
of the particle r 0 , the binary mass ratio M p /M c , and the binary eccentricity e.y,. In all cases 
the initial orbital eccentricity of the particle is taken to be zero, eo = 0. Each column of the 
figure corresponds to runs where we vary one of these parameters and fix the other three at 
fiducial (reference standard) values. The fiducial parameters that we adopt for this figure 
are io = 60°, r 0 = 0.2ab,eb = 0, and M p /M c = 1. The fiducial value of r 0 is chosen as a 
characteristic radius that typically lies within the disks that we simulate in Section [3j 


The left hand panels show that when the orbit of the particle is highly inclined 
(well above critical KL angle i cr ; see the green and blue lines), the KL cycle period is 
relatively insensitive to inclination, as expected by Equation (J3]) that has no dependence 
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Fig. 1.— Eccentricity (upper row) and inclination (bottom row) evolution of a test particle 
for different values of the initial orbital inclination xq, initial orbital radius tq, binary mass 
ratio M p /M c , and binary eccentricity e^- The binary orbit has a semi-major axis of a\ } . Time 
is in units of binary orbital period Py,. These results are obtained from three-body numerical 
simulations. (Color online) 























10 


on inclination. The simulated period is sometimes close to that predicted by Equation (J3|) . 
For i 0 = 70° (blue line), the eccentricity oscillation has e max ~ 0.93 and tkl — 15Pb, while 
Equations (J3]) and (d]) predict that e max = 0.90 and r KL = 16Pb- When the initial orbital 
tilt is just above the critical KL angle, the KL cycle period given by Equation (J3]) becomes 
less accurate. For io = 45° (red line), the period from the simulation is tkl — 33Pb, which 
is more than twice the analytical value. 

The second column shows the effect of changing the initial orbital radius of the particle 
r 0 . The KL periods for r 0 = 0.15ab, 0.2ab, and 0.25ab are ~ 30Pb, ~ 15Pb, and ~ 9Pb, 
respectively. These values agree well (within ~ 25%) with the predictions of Equation (J3J) 
that scales as tkl oc 7% 3 ^ 2 and produces KL periods of 24Pb, 16Pb, and UPb, respectively. 
Note that with r 0 = 0.25ab, the KL cycle period becomes shorter and the oscillation 
amplitude is reduced at later time (e m i n does not reach zero). This is likely because at 
this orbital radius, the perturbing potential from the companion is more significant. As a 
result, the stronger influence makes the analytic description as a secular (gradual) effect 
less accurate and causes the oscillations to be irregular. 

The effect of varying the binary mass ratio is shown in the third column. Equation (13|) 
predicts that r K L/Pb oc ( M c /M p + 1 ) 1//2 (M C /M p ) 1 ' /2 such that r K l = 10Pb, 16Pb, and 27Pb 
for binary mass ratios M p /M c = 2, 1, and 1/2, respectively. The results of the numerical 
simulations give us Tkl = 7Pb, 15Pb, and 32Pb. The agreement is then within about 30%. 

In the fourth column, we show the effect of varying the binary eccentricity. As 
discussed above, previous studies have shown that for eccentric orbit binaries, the maximum 
eccentricity can grow beyond e max given by Equation (J4]) and the maximum inclination can 
grow above the initial value of i 0 . As expected, for a somewhat higher binary eccentricity 
eb = 0.2 (green lines), the maximum and minimum particle eccentricity sometimes increases 
somewhat in each KL oscillation. The oscillation period is shorter than in the circular case, 
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as predicted by Equation (f3j), since the companion object interacts more strongly with the 
particle at binary periastron. According to Equation d3J); for eb = 0.2, the KL oscillation 
period is approximately 14Pb ; whereas the numerically determined period in Figure |T] is 
about 10Pb- Therefore, the reduction in the KL oscillation period due to binary orbital 
eccentricity is somewhat more severe than what is predicted analytically . 


3. KL cycle of a hydrodynamic disk 


use the SPH code PHANTOM 

(Loda 

ito & Price 

2010; 

Price & Federrath 

2010; 

Nixon 

2012; 

Nixon & King 

2012; 

Price 

2012; 

Nixon et al 

2013 

). The disk is initially circular, centere 


around the binary component of mass M c , and subject to perturbations by the other binary 
component of mass M p . The disk orbital plane is initially misaligned with respect to the 
binary orbital plane. The disk initially extends from radius r in to radius r out . The inner 
boundary of the simulated region is set to the disk initial inner disk radius r in . As particles 
move to r < r; n , they are removed from the simulation. We also impose an inner boundary 
radius around the perturbing companion, since some disk mass can be transferred to that 
component. The PHANTOM code adopts cubic spline kernel as the smoothing kernel. The 
softening length in the gravitational force is taken to be the same as the smoothing length. 
The number of neighbors is roughly constant at N neig h ~ 58. 

We define a certain set of fiducial (reference) model parameters. The disk is locally 
isothermal with sound speed c s oc r -3 / 4 and disk aspect ratio H/r = 0.035 at the initial 


inner disk radius r- m . These parameters all ow both a and t 


be constant throughout the disk radius flLodato fe Pringle 


le smoothing length ( h) /H to 


200 7). We employ an explicit 


accretion disk viscosity wh ich corresponds to an approximate 


Shakura & Sunvaev 

1973 

) over the disk (see 

Lodato & Price 

2010) 


y con stant a parameter 


2010). The viscous stresses 
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include a nonlinear term with a coefficient (3av — 2 (AV stands for artificial viscosity) that 
suppresses interparticle penetration, as is standard in SPH codes. The fiducial binary mass 
ratio is unity and the binary orbital eccentricity is zero (circular orbit). This fiducial setup 
is similar to that presented in Paper I. The only difference is that the initial number of 
particles in the fiducial model is set to N = 3 x 10 5 , while Paper 1 used N = 1 x 10 6 . 
The lower value adopted here reduces computational overhead in carrying out the various 
simulations, while providing results that are well converged in N, as we discuss below. 


Table [T] lists the simulation parameters. We consider several models whose parameters 
deviate from the fiducial model, but not all combinations of parameters were simulated. 
Instead, we consider variations of one parameter with all others fixed at the fiducial values, 
as we did for particle orbits in Figure [0 For different mass ratios, we also adjust the initial 
inner and outer disk radii. We scale the initial disk inner radius from the central binary 
component by a factor of (M C /(M C + Mp)) 1 / 3 . As discussed above, the accretion radius 
of the central object is taken to be the same as the initial inner disk radius. The inner 
boundary condition is that any SPH particle that goes inside accretion radius is removed 
from the simulation. The inner boundary radius of the perturbing component is set to 
r in (M p /M c y/ 3 . The initial disk outer radius r out is chosen to be that of a tidally truncated 


disk in a coplanar binary (jPaczyriski 


1977). Thus, when we change the mass ratio of the 


binary, we also change the initial disk outer radius. However, in a misaligned binary the 
tidal torque on the disk is weak ened and the disk can expand somewhat beyond this radius 


( iLubow. Mart i n fc Nix on 


2015). We ignore the effects of disk self-gravity in the simulations 


presented here. That means the disk does not feel its own gravity, but the stars can feel it. 
We set the ratio of the initial disk mass to the binary total mass to a very low value of 0.1% 
such that the binary orbit is hardly affected by disk’s gravity. This mass corresponds to a 
large initial Toomre Q of ~ 30. We will devote a separate paper to studying the effects of 
higher disk mass and stronger self-gravity. 








Table 1. Parameters of the SPH simulations for Binary Systems with Semi-major Axis of 




Binary and Disk parameters 

Symbol 

Fiducial Value 

Values 

Mass ratio of binary components 

Mp/Mc 

1 

[0.25, 0.5, 1, 2] 

Binary orbital eccentricity 

eb 

0 

[0, 0.2, 0.5] 

Initial number of particles 

7V/10 5 

3 

[2, 3, 5, 10] 

Initial disk mass 

M disk /(M C + Af p ) 

0.001 

0.001 

Initial disk outer radius 

^out/ 

0.25 

[0.15, 0.25] 

Initial disk inner radius 

Tin/ab 

0.025 

[0.015, 0.02, 0.025] 

Mass accretion radius 

^acc/&b 

0.025 

[0.015, 0.02, 0.025] 

Disk viscosity parameter 

a 

0.1 

[0.01, 0.1] 

Disk aspect ratio 

H/r (r = r in ) 

0.035 

[0.02, 0.035, 0.05, 0.065] 

Initial disk surface density profile £ oc r _T 

7 

1.5 

[0.5, 1.0, 1.5] 

Initial disk inclination 

*0 

60° 

[45°, 50°, 55°,60°] 
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In post-processing the simulations, we compute the orbital eccentricity and inclination 
for each particle using its position and velocity information according to Equations (JTj) 
and (ED- We divide the disk into 100 radial bins and calculate the mean properties of the 
particles within each bin, such as the inclination and eccentricity. The properties of the 
particles in these radial bins are not fully independent of each other because the disk is 
generally eccentric. Particles residing one bin follow streamlines that cover other nearby 
bins. Therefore, some correlation of mean properties occurs in nearby bins. In most of 
our analysis, we apply fixed radial bins that are fairly widely separated to reduce these 
correlation effects. An alternative approach is to apply bins based on particle semi-major 
axis a. We have found that binning in r and a give similar results for most of our analysis 
and so we generally report results with bins in r for simplicity. For the determination of 
disk warping, however, we adopt bins based on particle semi-major axis, as described later. 


The KL period of a test particle (Equation ([3])) depends on the particle’s orbital radius. 
Following Paper I, we determine an analytic estimate of the KL global period of a disk by 
taking an angular momentum weighted average involving the local KL period of particle 
orbits as follows 



(9) 


L 


rout r KT 1 Er 3 A /^dr 


where Tkl in the denominator is given by Equation ()3|). This estimate is equivalent to 
determining the disk global precession frequency as the integrated torque divided by 
the total disk angular momentum. This estimate is crude because we apply the angular 
momentum of circular orbits, while the disk is generally eccentric. In addition, the analytic 
particle orbit period in Equation (J3]) is itself somewhat approximate, as we found in Section 


El 


Figure El shows the evolution of the disk eccentricity and inclination at three different 
radii in the disk for various initial numbers of SPH particles. In this figure and following 
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similar figures, disk quantities are taken from three individual radial bins spanning 0.1, 0.2 
and 0.3afc, respectively. Also, we plot results up to a time of t — 28Pb- At this time, the disk 
has typically lost 90% of its initial particles. Most of the lost particles are accreted onto the 
central object. A small fraction of the particles are ejected from the disk, typically less than 
10%. Most of ejected particles are accreted onto the companion and others remain in a 
circumbinary orbit. From Figure [21 we see that the effect of the numerical resolution on the 
disk KL mechanism is in general quite small, especially for the first two cycles. Although 
the cycle period at later times in the low resolution run is slightly shorter than that in high 
resolution run, all the curves damp to nearly the same level of eccentricity e ~ 0.2 and 
inclination i ~ 28° at a time t = 28Pb- Given this insensitivity of the overall disk KL cycle 
pattern to the numerical resolution, we adopt N = 3 x 10 5 as the fiducial value that we use 
throughout the rest of the paper. 

Figure [2] shows some differences from the test particle simulation results in Figure [D 
We consider how Equations ([I]) and ([2]) might apply to a disk. We may expect a and e 
to represent a measure of the disk’s global semi-major axis and eccentricity, respectively, 
at some representative radius in the disk, such as the r = 0.2ab that we applied in Figure 
[0 Locally, a disk can transport angular momentum through waves and viscous torques. 
However, globally the vertical (perpendicular to the binary orbit plane) component of the 
disk angular momentum is conserved in the presence of an axisymmetric ring, apart from 
losses due to particles leaving the disk. On the other hand, the disk can dissipate energy, 
unlike the case of test particles. Therefore, the evolution of the disk is expected to depart 
from that of test particles, as given by Equations ([I]) and (J2D- In Figure [21 we see that 
the KL oscillations in the fluid case undergo damping in which the maxima of eccentricity 
decrease in time. In addition, the eccentricity minima do not drop to zero as in the particle 
case. The inclination drops below the minimum tilt of 39° for test particle oscillations. 
After the KL oscillations substantially damp, there is a significant residual disk eccentricity 


16 


e ~ 0.2. 


""he disk viscous timescale is estimated from Equation (18') of 


Lvnden-Be ll & Pringle 


(1197411 as t v rs./ r 2 /(120 20Pb for the fiducial model at r = 0.25ab. We have also verified 


that the disk density evolution of the fiducial model roughly foll ows t he results obt ained 


from of a ID viscous disk model (omits pressure effects) (e.g., 


Lvnden-Be ll & Pringle 


1974 ) 


for the same initial disk profile and level of viscosity. Therefore, over the course of the 
simulations, we expect the density profile to evolve from the initial one. 


The first peak value of the eccentricity seen in the upper row of the plot e ~ 0.70 is 
approximately consistent with the test particle results shown in Figure Q] for maximum 
e ~ 0.78 with an initial tilt of 60°, although the test particles have a somewhat shorter KL 
period. 

Figure [3] shows the effects of the disk viscosity on the KL mechanism. For lower 
viscosity (red curve; a = 0.01), the time at which the first eccentricity cycle peaks is 
basically unaffected. However, both the eccentricity and the inclination damp at a lower 
rate because the disk is less dissipative. For a = 0.01 at later times, the midpoints of 
the inclination KL oscillations settle to a tilt that is near the critical value of 39° for KL 
particle oscillations. Note that there is a substantial remnant disk eccentricity e ~ 0.3 after 
the KL oscillations significantly damp. The bottom panel of the figure plots the expected 
KL oscillation periods based on Equation ([9]). The estimated period is ~ 15Pb that is to be 
compared with the period obtained in the SPH simulations (upper panel) of ~ 7Pb- The 
agreement is crude (factor of 2), likely for reasons noted below Equation ©. 

The evolution of the disk eccentricity and inclination are well coordinated across 
different radii of the disk, as is evident from the similarity of the eccentricity and inclination 
curves across different radii in Figures [2] and [3l Figure [4] shows an edge-on view of the 
fiducial disk at t = 0 (left panel) and at t = 13Pb (right panel) when the line of nodes 
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r=0.1a 


r=n 9a 


r=H 9a 



10 15 

t/P b 


10 15 

t/P b 


Fig. 2.— Eccentricity (upper row) and inclination (bottom row) evolution of the disk at 
three different radii from the central star for different initial numbers of particles, db is the 
binary semimajor axis and Pb is the binary orbital period. The other parameter values are 
given in the third column of Table [Q (color online) 
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Fig. 3.— The upper figure is similar to Figure [2j but showing the effect of varying the disk 
viscosity parameter a. The lower figure shows the estimated disk KL period as a function 
of time (according to Equation (l9jl ). (color online) 
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Fig. 4.— View of the fiducial disk toward the x — z plane at times of 0 (left) and 13 Pb 
(right). The binary orbit is in the x — y plane (i.e., the perturbing object moves into and 
out of the page). The central mass is denoted by the black dot. The color coding is for the 
logarithm base 10 of the column density (i.e., density integrated along the line of sight) in 
units of (M c + M p )/a^. (Color online) 
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Fig. 5.— Left panel: the magnitude of the difference between the average disk angular 
momentum unit vectors for particles that lie in narrow bins of semi-major axis values centered 
on a = O.lcib and a = 0.2ctb, averaged in time over one binary orbit, plotted as a function of 
time for runs in Figure El Right panel: the difference between average disk tilts for particles 
in narrow bins of semi-major axis values centered on a = O.lab and a = 0.2etb, averaged in 
time over one binary orbit, plotted as a function of time for runs in Figure |3l 
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Fig. 6.— Similar to Figure |3l but showing the effect of varying the initial disk aspect ratio, 
(color online) 
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Fig. 7.— Similar to Figure El but showing the effect of varying the initial disk surface 
density profile. The lower right figure plots the azimuthally averaged surface density profiles 
at t — OPb and t = llPb- The surface density is in units of ( M c + M p )/a^ and r is in units 
of ab- (color online) 
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Fig. 8.— Similar to Figure [3l but showing the effect of varying the initial disk inner radius, 
(color online) 
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Fig. 9.— Similar to Figured but showing the effect of varying the initial disk outer radius, 
(color online) 
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has rotated by about 180°. This figure shows that the disk remains fairly flat, but quite 
lopsided, after substantial nodal and KL phase oscillation changes. 

We consider how inclination is radially communicated across the a disk. A KL disk 
undergoes tilt oscillations, along with nodal precession. A disk that is in good radial 
communication would be expected to remain nearly flat, that is, have a small degree of 
warping. The disk warping is often measured by rd r £ where £ is a unit vector that lies 
perpendicular to the local disk orbital plane. In the application to our simulation results, 
£(r) would be determined by the average value of £ for particles in a narrow radial bin 
centered on radius r. However, as noted above, there are artificial correlations in quantities 
such as £(r). The reason is that eccentric streamlines carry a single value of £ and extend 
radially in the disk. A better measure of warping is ad a £ where £(a) is determined by the 
average value of £ for particles whose semi-major axes lie in a narrow bin of size 3.5 x 10~ 3 ab 
centered about semi-major axis a. The reason for using a is that (noncrossing) streamlines 
in a fluid disk can be uniquely labeled by their a values. Therefore, ad a £ measures warping 
across streamlines. 

As an approximate measure of disk warping, we compute 

<5j = |j2-ji|, (10) 

where ji and j 2 are the average of the disk unit angular momentum vectors in bins of 
particle semi-major axis centered on aq and a 2 , respectively, that are separated by aq apart. 
We choose these values to be aq = O.lab and a 2 = 0.2ab- For the KL effect on test particles 
in a circular orbit binary, a is conserved (Equation (JT]) ). By a time ~ 10Pb, the value of Sj 
for test particles is of order unity, since the outer particle that begins on a circular orbit 
with r 0 = a = 0.2 a b undergoes a substantial change in tilt, as seen in Figure [I] as well as 
substantial nodal precession. 

Radial communication in a disk is diffusive for a > H/r and wave-like for a < H/r. 
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In the diffusive regime, the degree of warping is estimated by equating the binary nodal 


precessional torque to the horizontal (parallel to the disk 


involves the so-called v 2 viscosity flNataraian & Pringle 


component of viscous torque that 


19981 ). The u 2 effective viscosity 


owes its existence to the internal flows driven by p ressu re 


turbulent) viscosity, as described by 


gradients and limited by (e.g. 


Papaloizou & Pringle 


(1983). We roughly estimate 


then that the warping due to nodal precession is 


Sj ~ A sin (?) a — — 


w r 


n 


H 


(ii) 


where A is a coefficient of order unity, hi is the circular frequency at the disk radius 
r = 0.2ab, and u> p is the nodal precession frequency. We have estimated that A ~ 0.3 


b y so 


ving the ID viscous warp evolution equation given by Equation (10) of 


King et ah 


020131) for the conditions in the fiducial model that has oj p ~ 0.04flb and find that Sj ~ 0.1, 


which is roughly similar in value to the plotted results in left panel of Figure [5] for a = 0.1 
(blue curve). 


The lower viscosity a = 0.01 case plotted in left panel Figure [5] lies in the wave-like 
regime. The communication is a consequence of radial pressure forces. In this simulation 
the radial sound crossing timescale is substantially shorter than the KL oscillation or nodal 
precessional timescale. We have that r/(c s pKi,) ~ 0.1. As described in Paper 1, be caus e this 


ratio is small the c 


Lubow fc Ogilvie 


isk undergoes global rigid tilt oscillations (Larwood & Papaloizou 


1997|; 


200l|). In this case, the level of warping (plotted as the red curve on the 


left panel of Figure [5]) is similar to or weaker than in the viscous regime (blue curve). 


The warping Sj depends on the variation of both the disk tilt and direction of the line of 
nodes. It is affected by both nodal precession and KL disk tilt oscillations. We consider here 
the variations of only the disk tilt. The KL inclination oscillations can lead to variations in 
disk tilt across streamlines in the disk. These variations can be expected because the KL 
oscillation timescale varies with particle orbit period P and therefore with particle orbit 
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semi-major axis a. According to Equation (l3lh the KL oscillation period should differ by 
about a factor of ~ 3 in a change from a = O.lab to a = 0.2ab. Radial communication via 
viscous and pressure effects as described above play a role in maintaining a constant tilt. 
As seen in the right panel of Figure [5j the change in tilt between a = O.lab and a = 0.2ab is 
small. The difference in tilt angles between them is typically less than 3°. This change is 
similar to the disk aspect ratio H/r = 0.035 ~ 2.0°, while the disk oscillates over a range 
of tilt angles ~ 20°. In comparing the right panel of Figure [5] with Figure [H we see that 
tilt difference is phased with the KL cycle. When the disk inclination is minimum and 
eccentricity is maximum, the outer inclination i 2 (that has a = 0.2ab) is smaller than the 
inner inclination i\ (that has a = O.lab) . 

The effect of varying the disk aspect ratio is shown in Figure [6] Most noticeably, as the 
disk aspect ratio decreases (from the blue curve to the red curve), the first disk KL cycle 
gets delayed to a later time. Because the kinematic viscosity u varies with ( H/r ) 2 , lower 
values of H/r lead to a slower disk expansion and therefore a weaker perturbation from the 
companion object. Lower v values also lead to slower damping in both eccentricity and 
inclination. The bottom panel of the figure shows the expected KL oscillation periods based 
on Equation (J9]). As we saw in the case of Figure [3j the accuracy of this equation is only to 
about a factor of 2. Based on the bottom panel, the ratio of KL periods for H/r = 0.065 to 
H/r = 0.02 is about 1.3, while the ratio in the SPH simulations (upper panel) is about 1.5. 

We present in Figure [7| three runs using different power-law indices for the initial disk 
surface density profile. The differences in the results for these three cases is found to be 
quite small. The bottom left panel of this figure shows that the expected KL periods based 
on Equation (J9]) are about the same after a time of about 10Pb, as is in agreement with 
the SPH simulations in the upper panels. The bottom right panel shows that after initial 
adjustment, the density distributions are very similar. 
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In Figure [HI we show the effects of changing the disk inner radius. The first disk KL 
cycle does not seem to be affected by the location of the disk inner boundary. However, the 
cycle period does get slightly longer as the disk inner radius decreases. The bottom panel 
of this figure shows that the expected KL periods based on Equation (jUJ) are in qualitative 
agreement with this trend. This effect is likely be the result of having more material 
located further from the companion, closer to the central object, resulting in weaker global 
perturbations. 

Figure O shows the effects of changing the disk outer radius. With smaller initial 
disk outer radius, the disk initially responds more slowly since the perturbations from the 
companion object are weaker. Thus, the effect in this case is simply a time delay in the 
cycles. The bottom panel of this figure shows that the expected KL periods based on 
Equation (|9l) are in agreement with this description. The case with the smaller initial radius 
has an initial KL period that is much longer than for disk of large initial radius. By a time 
of ~ 15/j,. both cases have the nearly same KL periods, as seen in the SPH simulations 
results of the upper panels. 

Figure [10] shows the effects of changing the initial disk inclination. As z 0 varies from 
60° to 45°, the KL cycles occur at a progressively later time with a longer cycle period 
and smaller oscillation magnitude. This trend seems to be in line with the test particle 
results discussed in the previous section, where the KL oscillation period increased with 
decreasing inclination. Equation (JTJ) predicts eccentricity maxima of 0.76, 0.67, and 0.56 for 
initial tilts of 60°, 55°, and 50°, respectively, while Figure [TUI exhibits maximum e values 
during the first KL cycle of 0.7, 0.6, and 0.4. In the test particle case, for an initial tilt of 
io = 45°, the eccentricity can reach values of 0.4, based on Equation (JTJ), or 0.3, based on 
simulations (top-left panel of Figured]). However, for this initial inclination angle, the KL 
mechanism barely operates in the disk whereas it is still relatively strong in the test particle 
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runs (see the left column of Figure [T]). Thus, the critical initial tilt angle for the onset of 
KL oscillations in a hydrodynamical disk is somewhat higher than for a test particle (39°). 
On the other hand, in the first disk KL cycle, the inclination can drop below 39°, which 
does not occur in the particle case. 

Figure fill demonstrates the effects of varying the binary mass ratio. For M p /M c < 1, a 
smaller mass ratio leads to a later occurrence of the KL cycle and also a longer cycle period 
(from the green curve to the red curve), which again is in line with the test particle results. 
The bottom panel of this figure shows that the expected KL periods based on Equation (J9j) 
at a time of 25Pb are about 30Pb, 20Pb, 15Pb, and 7Pb for mass ratios of 1/4, 1/2, 1, and 2, 
respectively, whereas the SPH simulation results show 15Pb, UPb, 8Pb, and UPb- For both 
approaches for mass ratios less than or equal to unity, the KL period decreases with mass 
ratio, as expected. The period for the corresponding cases are again different by about a 
factor or 2. In the case where the perturbing companion is more massive with a mass ratio 
of 2, the period in the SPH simulation does not follow this trend and instead increases. 


The effects of varying the binary eccentricity are shown in Figure 033 The initial disk 
outer radius is r out = 0.2% instead of the fiducial value of 0.25%. This smaller initial disk 
outer radius was used because the increased binary eccentricity causes a stronger response 
in the disk and a smaller disk truncation radius, as is known to occur for the coplanar case 


(Artymowicz & Lubow 


3 ) 


1994]). Other parameters are set to their fiducial values. We do not 


plot the results for eb = 0.5 at r = 0.2% and r = 0.3% because the disk is truncated within 
this radius due to the smaller periastron distance. The KL oscillation periods are about the 
same in the three cases, within about 10% of each other. The maximum eccentricities are 
also about the same in the three cases, as we found for test particles in Section [2j In the 
eb = 0.5 case, the inclination decays to a value close to i CI = 39°, while the lower binary 
eccentricity cases decay to a smaller tilt angle. 
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4. Discussion and Summary 


In this paper, we investigated the conditions under which KL oscillations can operate 
on a fluid disk, as was first described in 


Martin et ah 


f)2014bl ). Such a disk undergoes 


a complicated evolution involving tilt and eccentricity oscillations, as well as the usual 
nodal and apsidal precession. We carried out this analysis by means of SPH simulations 
that include the effects of disk pressure and viscosity. We considered a range of disk 
viscosities (Figure E]), sound speeds (Figure E]), initial density profiles (Figures [71 [SI EJ), 
initial inclinations (Figure [TO]), binary mass ratios (Figure [TIT), and binary eccentricities 
(Figure [IS]). We found that the KL effect generally operates over a range of of disk and 
binary properties. 


The general picture we have is that an initially circular circumstellar disk that is 
sufficiently inclined (~ 45° - 135°) with respect to the orbit plane of a binary undergoes 
damped KL eccentricity and tilt oscillations. During these oscillations, the disk can develop 
substantial eccentricities, e > 0.5. Once the KL oscillations have mostly damped, the disk 


achieves a tilt in the prograde disk case (initial tilt less than 90°) i ~ 39° or somewhat 


that long, we expect that the disk eccentricity will eventually damp 
and the disk will become coplanar with the binary orbit plane (e.g., 
to viscous effects. 


(e.g., 

Ogilvie 

1001 

) 

King et al. 

201 

3) due 


The properties of KL disk oscillations are somewhat similar to the properties of KL test 
particle oscillations (Figure [1]). In both cases, inclination and eccentricity are interchanged. 
The peak eccentricity achieved by a disk in the first KL oscillation is generally close to 
the eccentricity achieved by test particles for the same initial inclination. In addition, for 
smaller mass perturbers, the KL oscillations are slower for both test particles and disks 
(Figure [111). 
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However, there are some significant differences between a fluid disk and test particles. 
The KL oscillation frequency of test particles varies with particle orbit semi-major axis 
as a 3 / 2 . On the other hand, a disk undergoes a large scale global oscillation with a small 
level of warping for typical parameters (Figures 0] and [5]). For a circular orbit binary, the 
inclination of prograde test particle orbits relative to the binary orbit plane must be greater 
than the critical angle i cr = arccos (a/3/5) ~ 39° for KL oscillations to occur. In the case of 
a disk, we find that the minimum inclination angle appears to be somewhat larger (> 45°) 
(see Figure flOl). For a circular orbit binary, the inclination of prograde test particle orbits 
oscillates between the initial value i 0 > i CI and i cr . For the case of a disk, however, the 
oscillation can dip below i cr (e.g., Figure [3]). Therefore, even though a disk requires a larger 
tilt to operate than a test particle orbit, the disk oscillations can reach lower tilt angles. 


Another major difference between test particles and disks is that disks can dissipate 
energy. As a result, the KL oscillations damp, but leave a residual disk eccentricity that 
we find to be typically a few tenths. We find that the damping rate depends on the 


level of disk viscosity (Figure EJ)- The disk eccentricities may a 


process involving a par ametric instability (IPapaloizou 


2005a 


so damp through a 


Og ilvi e fe Barker 


Barker & Ogilvie 


decay 


2014 ; 


2014). Our simulations may not be capable of resolving this instability 


that occurs on scales smaller than the disk thickness. It is also possible that the initial KL 
oscillation damping involves shocks, since the streamlines may attempt to cross at high 
eccentricity. We have found indirect evidence for such effects in some simulations that we 
carried out with a low value of viscosity parameter /?av- As mentioned in Section |3l this 
parameter is used in the SPH code to prevent interparticle penetration that artificially 
occurs in the code where shocks are present. For small values of (3av, the SPH simulations 
displayed unphysical behavior, symptomatic of interparticle penetration. Over longer 


timescales in the post-KL oscillation phase, the dis 
the binary orbit plane through viscous effects (e.g., 


£ is expectec 


King et al. 


to eventually align with 


2013) 
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For an eccentric orbit binary, the KL test particle oscillations can reach to higher 
inclination angles than the initial tilt. Inclination an gles can reach 90° and even higher, 


resu 


t ing flip pi ng from prograde to retrograde orbits (ILithwick fe Naoz 


2011 


2013a 


b; 


Liu et ah 


2011 


Naoz et ah 


2015 ). In reaching 90°, the orbit becomes radial. KL disk 


oscillations cannot achieve a similar increased level of inclination (certainly not to 
eccentricity close to unity) because of dissipation. 


In this paper, we have omitted the effects of dis k se 


KL oscillations for sufficiently massive disks ([Batygin 


f-gravity. Such effects may 


2012 


Martin et ah 


suppress 


2014b). We 


intend to describe the influence of self-gravity in a future paper. 

Confirmation of the KL mechanism acting on a disk would come from an observation 
of an eccentric and misaligned disk in a binary system. The KL oscillations are applicable 
to binary systems on all scales. For example, the oscillations could play a role in planet 
formation and evolution around one component of a binary star. There is observational 


evidence^ th at s ome exoplanets may be unc 


Wu & Murray 

2003; 

Takeda & Rasio 

2005) 


ergoing KL cycles in wide binaries (e.g. 


2005 ). In future work, we plan to investigate the 


evolution of a planet-disk system that orbits around one component of a binary system in 
order to understand how these planets can form. 


Be/X-ray 


(Martin et al. 


oinaries may also be influenced by the KL oscillations of a misaligned disk 


2014aUbl) . Be stars m ay form disks from material that is e 


equator of the rapidly rotating star (ILee. Osaki fc Saio 


1991 


Hanuschik 


ectec 


1996 


from the 


Carciofi 


2011). Be/X-ray binaries contain a Be star with a rotation axis that is misaligned to that 


of the orbit of a companion neutron star (see Table 1 in 


Martin et al. 


2011 


, for some 


observations). We found in our simulations that roughly 10% of the initial disk mass 
is captured by the companion object for the fiducial model. As a disk undergoes KL 
oscillations to maximum eccentricity, its apastron radius grows which allows this capture to 
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occur. This mass exchange may explain the so-called giant 


the accretion of disk mass by the companion neutron star (Martin et al 


Type II outbu rsts as being due 


2014 


aj)- 


Finally, the KL oscillations of a disk may also be relevant to SMBH binaries. It 
i s possib le that the circumbinary disk around a SMBH binary will be misaligned (e.g. 


Nixon ct al. 


2011a.b). Thus, as gas accretes inwards from the circumbinary disk, misaligned 


disks will form around each SMBH. These disks will be unstable to KL oscillations if 
the misalignment angle is sufficiently large. This effect could have implications for the 
subsequent disk evolution and star formation in these systems. 
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2007) for data analysis and rendering of figures. 
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Fig. 10.— Similar to Figured but showing the effect of varying the initial disk inclination 
angle, (color online) 
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Fig. 11.— Similar to Figure [3l but showing the effect of varying the binary mass ratio. Note 
that here we show results up to t — 40Pb in the upper figure, (color online) 
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Fig. 12.— Similar to Figure El but showing the effect of varying the binary eccentricity, 
(color online) 
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